What do we mean by a “pipeline” or “workflow”?

What do we mean by a “pipeline” or “workflow”?

We may colloquially refer to any consecutive series of steps in an analysis as a pipeline or workflow.

But here, I am using these words to mean something that can be executed from start to finish with a single command.

A basic pipeline script

samples=(smpA smpC smpG)

# Trim:
for sample in ${samples[@]}; do
    scripts/trim.sh data/"$sample".fq > res/"$sample"_trim.fq
done

# Map:
for sample in ${samples[@]}; do
    scripts/map.sh res/"$sample"_trim.fq > res/"$sample".bam
done

# Count:
scripts/count.sh ${samples[@]} > res/count_table.txt

Why create a pipeline?

What are the advantages of creating such a pipeline, rather than running scripts one-by-one as needed?

• Rerunning everything is much easier.

• Re-applying the same set of analyses in a different project is much easier.

Challenges with basic pipeline scripts

How to rerun parts of the pipeline flexibly? This may be necessary after, e.g.:

• Some scripts fail for all or some samples.

• Adding a sample.

• Needing to modify a script or settings somewhere halfway the pipeline.

for sample in ${samples[@]}; do
    sbatch scripts/trim.sh data/"$sample".fq > res/"$sample"_trim.fq
    sbatch scripts/map.sh res/"$sample"_trim.fq > res/"$sample".bam
done

sbatch scripts/count.sh ${samples[@]} > res/count_table.txt

Potential solutions to these challenges

• Informal pipeline
  Rather than running the pipeline script as an actual script, just use it as scaffolding and submit interactively jobs one-by-one / group-by-group. You don’t have a true pipeline but that may be OK.

• Push the limits of the Bash and Slurm tool set
  Work with if statements, many arguments to scripts, and Slurm “job dependencies” — but this is very hard to manage for more complex workflows.

• Use a formal workflow management tool.

The need for specialized tools

Perkel 2019 - https://www.nature.com/articles/d41586-019-02619-z

Workflow management systems

Pipeline/workflow tools, often called “workflow management systems”, provide ways to formally describe and execute pipelines.

Some of the most commonly used (command-line based) options in bioinformatics are:

• Nextflow
• Snakemake
• Common Workflow Language (CWL)
• Workflow Description Language (WDL)

Advantages of formal workflow management

Automation

• Detect & rerun upon changes in input files and failed steps.
• Easily run for other data sets.
• Automate Slurm job submissions.
• Integration with software management.

Writing pipelines vs. using them

Most workflow tools are small “domain-specific” languages -DSLs-, often a sort of extension of a more general language (Python for Snakemake, Groovy for Nextflow).

nf-core curates Nextflow pipelines

The “nf-core” initiative (https://nf-co.re) curates a set of best-practice, flexible, and well-documented pipelines written with Nextflow.

It also has some of its own tooling built on top of Nextflow to help contributors create more robust and user-friendly pipelines.

The nf-core paper

https://www.nature.com/articles/s41587-020-0439-x

The Nextflow paper

https://www.nature.com/articles/nbt.3820

nf-core pipelines

nf-core currently has 58 complete pipelines — these are the four most popular ones:

Today’s lab: nf-core rnaseq

We will run the nf-core rnaseq pipeline (https://nf-co.re/rnaseq), a reference-based RNA-seq pipeline.

Notes on running an nf-core pipeline

Running a pipeline like this is a bit different (and more involved) than running a typical piece of bioinformatics software — it is, after all, stitching together many operations. Some considerations:

• We only need to install Nextflow and download the workflow files.

• The pipeline runs all its constituent tools via Singularity containers (this is most common, and recommended) or via Conda environments.

RNA-seq to do what?

The nf-core rnaseq pipeline is meant for RNA-seq projects that:

• Attempt to sequence only mRNA while avoiding non-coding RNAs (“mRNA-seq”).
• Do not distinguish between RNA from different cell types (“bulk RNA-seq”).
• Use short reads (≤150 bp) that do not cover full transcripts but do uniquely ID genes.
• Use reference genomes (are reference-based) to associate reads with genes.

Two main parts to RNA-seq data analysis

There are two main parts to the kind of RNA-seq data analysis I just described:

• From reads to counts
  
  • Generating a count table using the reads & the reference genome.
  • This what nf-core rnaseq does, which we will run in today’s lab.

• Count table analysis
  
  • Exploratory data & differential expression analysis is covered in the self-study lab.
  • Functional enrichment (GO, KEGG) analysis is typically also performed.

Overview of the pipeline again

The inputs are FASTQ and reference genome files (assembly & annotation), and the outputs include a gene count table and many “QC outputs”.

The main steps

Read QC and pre-processing

1. Read QC (FastQC)

2. Adapter and quality trimming (TrimGalore)

3. Optional removal of rRNA (SortMeRNA) — off by default, but we will include this

Alignment & quantification

4. Alignment to the reference genome/transcriptome (STAR)

5. Gene expression quantification (Salmon)

Post-processing, QC, and reporting

6. Post-processing alignments: sort, index, mark duplicates (samtools, Picard)

7. Alignment/count QC (RSeQC, Qualimap, dupRadar, Preseq, DESeq2)

8. Create a QC/metrics report (MultiQC)

Pipeline variants

This pipeline is quite flexible and you can turn several steps off, add optional steps, and change individual options for most tools that the pipeline runs.

• Optional removal of contaminants (BBSplit)
  Map to 1 or more additional genomes whose sequences may be present as contamination, and remove reads that map better to contaminant genomes.

From reads to counts: overview

This part is “bioinformatics-heavy” with large files, a need for lots of computing power, and command-line (Unix shell) programs — it specifically involves:

1. Read pre-processing

2. Aligning reads to a reference genome

3. Quality control (QC) of both the reads and the alignments

4. Quantifying expression levels

Reads to counts: read pre-processing

Read pre-processing includes the following steps:

• Checking the quantity and quality of your reads
  • Does not change your data, but helps decide next steps / sample exclusion
  • Also useful to check for contamination, library complexity, and adapter content

Reads to counts: alignment to a reference genome

The alignment of reads to a reference genome needs to be “splice-aware”.

Berge et al. 2019

Reads to counts: alignment to a reference genome

Alternatively, you can align to the transcriptome (i.e., all mature transcripts):

Berge et al. 2019

Reads to counts: alignment QC

Here are some examples of the kinds of things to pay attention to:

• Alignment rates
  What percentage of reads was successfully aligned? (Should be >80%)

Reads to counts: quantification

At heart, a simple counting exercise once you have the alignments in hand.
But made more complicated by sequencing biases and multi-mapping reads.

Count table analysis

The second part of RNA-seq data analysis involves analyzing the count table.
In contrast to the first part, this can be done on a laptop and instead is heavier on statistics, data visualization and biological interpretation.

It is typically done with the R language, and common steps include:

• Principal Component Analysis (PCA)
  Assessing overall sample clustering patterns

• Differential Expression (DE) analysis
  Finding genes that differ in expression level between sample groups (DEGs)

• Functional enrichment analysis
  See whether certain gene function groupings are over-represented among DEGs

Minimal Nextflow pipeline example

params.greeting = 'Hello world!' 
greeting_ch = Channel.of(params.greeting) 

process SPLITLETTERS { 
    input: 
    val x 

    output: 
    path 'chunk_*' 

    script: 
    """
    printf '$x' | split -b 6 - chunk_
    """
} 

process CONVERTTOUPPER { 
    input: 
    path y 

    output: 
    stdout 

    script: 
    """
    cat $y | tr '[a-z]' '[A-Z]'
    """
} 

workflow { 
    letters_ch = SPLITLETTERS(greeting_ch) 
    results_ch = CONVERTTOUPPER(letters_ch.flatten()) 
    results_ch.view { it } 
}